Limited Ferromagnetic Interactions in Monolayers of MPS3 (M = Mn and Ni)

We present a systematic study of the electronic and magnetic properties of two-dimensional ordered alloys, consisting of two representative hosts (MnPS3 and NiPS3) of transition metal phosphorus trichalcogenides doped with 3d elements. For both hosts, our DFT + U calculations are able to qualitatively reproduce the ratios and signs of all experimentally observed magnetic couplings. The relative strength of all antiferromagnetic exchange couplings, both in MnPS3 and in NiPS3, can successfully be explained using an effective direct exchange model: it reveals that the third-neighbor exchange dominates in NiPS3 due to the filling of the t2g subshell, whereas for MnPS3, the first-neighbor exchange prevails, owing to the presence of the t2g magnetism. On the other hand, the nearest neighbor ferromagnetic coupling in NiPS3 can only be explained using a more complex superexchange model and is (also) largely triggered by the absence of the t2g magnetism. For the doped systems, the DFT + U calculations revealed that magnetic impurities do not affect the magnetic ordering observed in the pure phases, and thus, in general in these systems, ferromagnetism may not be easily induced by such a kind of elemental doping. However, unlike for the hosts, the first and second (dopant–host) exchange couplings are of similar order of magnitude. This leads to frustration in the case of antiferromagnetic coupling and may be one of the reasons of the observed lower magnetic ordering temperature of the doped systems.


■ INTRODUCTION
Nonmagnetic van der Waals layered materials such as transition metal dichalcogenides have been extensively studied over the last several years. 1 Just recently, the intrinsic ferromagnetism in the true 2D limit has been reported, 2 initiating increasing excitement in the spintronics and 2D material communities. In particular, the long-range magnetic order has been observed for insulating CrI 3 , 3 semiconducting Cr 2 Ge 2 Te 6 , 4 and metallic Fe 3 GeTe 2 5 compounds. In addition, topological spin structures have been predicted for 2D materials. 6 Currently, the attention is on transition metal phosphorus trichalcogenides (MPX 3 , where M stands for a transition atom and X = S and Se), which could be easily exfoliated down to monolayers 7,8 and are semiconducting materials with a wide range of band gaps. 9 The MPX 3 structures exhibit various antiferromagnetic (AFM) arrangements within the magnetic ions, which are theoretically expected to be measurable using different light polarizations. 10 Interestingly, the metal to insulator transition and superconductivity phase have been observed in this compound family. 11 −13 In particular, the AFM insulator phase in bulk FePS 3 can be melted and transformed into the superconducting phase under high pressure, providing similarity to the high-T c cuprate phase diagram. 14 In addition, theoretical predictions point to the existence of a large binding energy of excitons in MnPS 3 , whereas the experimental reports have observed excitons in few layers of NiPS 3 strongly related to magnetic order. 13,15,16 Recent reports have demonstrated an all-optical control of the magnetic anisotropy in NiPS 3 by tuning the photon energy in resonance with an orbital transition between crystal field-split levels. 17 The aforementioned demonstrates that this family of compounds is an ideal platform to study correlation effects in the true 2D limit.
In contrast to ferromagnetic (FM) materials, the antiferromagnets exhibit limited applications, mostly in the terahertz regime as ultra-fast components or specialized embedded memory−logic devices. 18−20 Most of the current applications of magnetic crystals are based on the FM semiconductors, for which the band gap and FM order are crucial factors. The magnetic-phase transitions for the MPX 3 materials can be accomplished by applying stress, 21 changing the carrier concentration or applying voltage. 22 In addition, the "M" atoms in MPX 3 crystals might be substituted with other transition metal atoms inducing the FM order, as recently reported for a particular concentration of the non-magnetic dopants in the CrPSe 3 host, resulting in a half-metallic FM state. 23 Moreover, a series of mixed systems in a bulk form 24−35 has been experimentally realized, thus entering a new playground of magnetic phases.
In addition, the magnetic and electronic properties of MPX 3 materials and critical Neél temperature (from T N = 78 K for MnPS 3 36−38 up to T N = 155 K for NiPS 3 36,37 ) strongly depend on the type of magnetic ion in the host. Thus, an elemental substitution could be an efficient way to tune the magnetic properties of atomically thin layers, by changing the lattice parameters and magnetic moments. These quantities result in manipulation of the exchange interactions that could be an effective way to engineer highly functional materials, similar to magnetic heterostructures. 39 In particular, the main reason behind the idea that adding ions with the partially filled dshells into the system might lead to the enhancement of the FM interactions is related to the so-called double-exchange mechanism: 40 the hopping between two correlated ions with different valences (and thus comparable Coulomb interactions) and relatively strong Hund's exchange is energetically favored provided that the spins on the neighboring ions are aligned in a parallel fashion. Such a mechanism might for instance be at play along the Cr−Mn bond: here the hopping from the manganese e g orbitals (e g 2 configuration) to the chromium e g orbitals (e g 1 configuration) lowers the total energy of the system provided that the spins on Cr and Mn are aligned ferromagnetically (for the antiparallel configuration, the strong Hund's exchange makes such a hopping energetically unfavorable). Note that such a double-exchange mechanism is similar to the one observed in the doped manganites which are FMalbeit here it is a bit more complex, for it involves two different magnetic ions (nevertheless, the latter should not be that important, for Cr and Mn have comparable values of the Coulomb interactions, cf. Table II of ref 41).
Altogether, in this work, we study the magnetic properties of the MPX 3 monolayers and try to understand (i) how FM interactions can be stabilized in these compounds and (ii) whether one can easily modify these compounds by elemental doping so that the FM interactions can be strengthened. To this end, we examine two representative hosts MnPS 3 and NiPS 3 in the undoped case and at a particular doping concentration of magnetic ions. Using first-principles calculations, we give qualitative and quantitative explanations of the origin of the exchange coupling strengths up to the third nearest neighbors and their respective ratios for MnPS 3 and NiPS 3 structures. Considering the model Hamiltonians with ab initio parametrization, we discuss the competition between the direct exchange and the superexchange mechanisms for the host structures. Next, we study various substantial sites of dopant atoms with mixed spins and mixed nearest-neighbor magnetic interactions. Here, again using first-principles calculations, we examine in detail the mixed exchange coupling parameters between the metal host and dopant atoms.
In the next section, we describe the computational details. In the third section, we will present our Results and Discussion while the last section is devoted to the Conclusions.

■ METHODS
DFT Computational Details. The first-principles calculations are performed in the framework of spin-polarized density functional theory (DFT) as is implemented in the VASP code. 42,43 The electron−ion interaction is modeled by using PAW pseudopotentials 44,45 with 3s and 3p states for P and S atoms and 3d and 4s states for Mn, Ni, and Cr being treated as valence states. The Perdew−Burke−Ernzerhof (PBE) exchange−correlation functional is employed. 46 The kinetic energy cutoff for the plane-wave expansion of the wave functions is set to 400 eV. A k-mesh of 10 × 6 × 2 is taken to sample an irreducible first Brillouin zone of the rectangular planar cell (see Figure S1 in Supporting Information) containing 20 atoms including four transition metal atoms. The lattice parameters have been fully optimized within the PBE + U approach for the magnetic ground state of the monolayers, assuming the rectangular supercell. In particular, the magnetic ground state (Hubbard U parameter) for MnPS 3 and NiPS 3 is AFM-N (U = 5 eV) and AFM-z (U = 6 eV), respectively. In the case of (2 × 2) supercell, which consists of four primitive hexagonal unit cells, the 5 × 5 × 1 k-mesh is chosen to obtain the optimized position of the atoms. Considering density of states (DOS) calculations, the denser k-mesh equal to 10 × 10 × 2 k-points is taken into account. The convergence criteria for the energy and force are set to 10 −5 eV and 10 −3 eV/Å, respectively. In order to properly model a monolayer system, 20 Å of a vacuum is added to neglect the spurious interaction between the image cells. Note that the standard exchange−correlation functionals are insufficient to account for a non-local nature of dispersive forces, which are crucial for layered materials and adsorption molecules on the surfaces. 47−50 Thus, the semi-empirical Grimme method 51 with a D3 parametrization is applied. 52 For the 2D materials, we can use the HSE for a better description of the gap; 53 however, in the case of magnetism also, GGA + U provides the same effect. We employ the GGA + U formalism proposed by Dudarev 54 to properly account for on-site Coulomb repulsion between 3d electrons of transition metal ions, by using effective Hubbard U parameters.
Note that the proper choice of the U values is not straightforward due to the lack of accurate experimental information on electronic properties. Also, the common choice to compare the band gaps obtained in DFT with experiments to judge the U values also need caution due to the fact that one-particle Kohn−Sham DOS cannot be directly compared to the measured data. Thus, we decided to compute the Hubbard U using the linear response method proposed by Cococcioni 55 for the monolayers of MnPS 3 and NiPS 3 , and we obtained 5.6 eV for the Ni and 5.3 eV for the Mn. Therefore, we have used U Cr = 4 eV, U Mn = 5 eV, and U Ni = 6 eV similar to the linear response results, which are typical U values reported for MPX 3 materials in previous reports. 21 Moreover, our Coulombic repulsions are close to the typical values of U in semiconductor compounds. Indeed, we find U = 6.4 for Ni 2+ and U = 3.5 eV for Cr in oxides, 56 while the typical value of U for Mn 2+ is U = 5 eV. 57,58 To calculate the AFM direct exchange for MnPS 3 and NiPS 3 , we extracted the real-space tight-binding Hamiltonian with atom-centered Wannier functions with d-like orbital projections on the transition metals using the Wannier90 code. 59,60 We calculated separately the hopping parameters for the orbitals, symmetric and antisymmetric with respect to the basal plane. The different symmetry and the separation in energy help to disentangle the two subsectors of the d manifold. 61,62 The calculation of the hopping parameters was carried out in the non-magnetic case to get rid of the magnetic effects and evaluate just the bare-band structure hopping parameters, and then, the hopping parameters will be used for the model Hamiltonian part. 63 In order to have parameters to use for the

■ RESULTS AND DISCUSSION
The results are presented as following: first, we present the results for the magnetic ground state of the hosts (pure MnPS 3 and NiPS 3 systems). In particular, we consider the exchange couplings within the DFT + U approach. Next, the AFM exchange mechanism is discussed within the minimal manybody model, and the AFM exchange coupling strengths are evaluated numerically using the Wannier basis with ab initio parametrization. Next, the qualitative explanation of the FM superexchange is presented. Finally, we present comprehensive studies of electronic and magnetic properties of benchmark alloys with a fixed concentration of dopants. The elemental substitution is employed at various atomic sites of the honeycomb lattice. Undoped Hosts MnPS 3 and NiPS 3 . Magnetic Couplings within the DFT + U Approach. The magnetic ground states of MnPS 3 and NiPS 3 exhibit AFM Neél (AFM-N)-and AFM zigzag (AFM-z)-type ordering (see Figure 1), respectively. The neutron diffraction data predicted that the Mn 2+ (3d 5 , S = 5/2) spins are slightly tilted (around 8°) from the perpendicular direction of the honeycomb lattice, 68 whereas the Ni 2+ (3d 8 , S = 1) spins are aligned within the honeycomb plane. In addition, due to the different filling of the 3d orbitals for various metals ions (Fe, Mn, Ni etc.), the size of the magnetic exchange coupling (J) also changes. The existence of magnetic ordering at finite temperature in 2D limit requires the magnetic anisotropy in accordance to the Mermin−Wagner theorem. Currently, only FePS 3 , which possesses a strong out-of-plane easy axis, has been experimentally reported to exhibit AFM order in the monolayer up to 118 K. 8 In order to explain the anisotropic order, the dipolar interactions between the magnetic moments or a single-ion anisotropy resulting from a nonzero spin−orbit coupling should be considered. The recent experimental reports, such as magnon band measurements, 69 support the claim that the dipolar interactions should be the leading term, whereas the electron spin resonance 70 and critical behavior measurements 71 indicate that the single-ion anisotropy might come into play. In addition, experimental observations demonstrated that the AFM ordering persists down to bilayer samples and is suppressed in the monolayer. 72 Notably, the recent theoretical report 73 has questioned the Raman criterion used for the monolayer studies therein, suggesting that NiPS 3 magnetic ordering could be presented in monolayer samples, as also indicated by strong two-magnon continuum existing in thin samples of NiPS 3 . 72 In the case of MnPS 3 , the magnetic order has been presented down to the bilayer and was reported to be absent in the monolayer. 74 The suppression of the Neél temperature in thin samples can be associated with reducing the interplanar coupling in atomically thin samples. 3,72 Aforementioned results demonstrate that the magnetic ordering in monolayers of MPX 3 is still under a hot debate and many experiments are being carried out to verify the theoretical predictions. Similar effects of strong spin fluctuation and absence of interlayer exchange coupling that weaken the long-range spin order in the 2D limit have been reported in other layered magnets such as Cr 2 Ge 2 Te 6 and CrI 3 . 3 Here, we focus on the rationally invariant Heisenberg contribution, whereas the dipolar and single-ion contributions are out of the scope of the present work. The latter has been estimated to 0.3 and 0.009 meV for NiPS 3 and MnPS 3 , respectively, and discussed theoretically in ref 66. Here, we present the results of exchange couplings up to the third nearest neighbor (for clarification, see Figure S1 in Supporting Information).
The exchange interaction up to the third nearest neighbor (J i M ) between the metal atoms of the hosts has been widely studied in a series of previous studies 21,66,67,75−77 (see Table  1). Note that the prediction of the magnetic ground state within the DFT calculations does not depend on Hubbard interaction, whereas it is well known that the exchange coupling strength is sensitive to both U and the lattice parameters. We set the Hubbard U parameter to U = 5 eV and U = 6 eV for the 3d orbitals of Mn and Ni atoms, respectively. These values are calculated from first principles using the Cococcioni approach. 55 Our predicted J i M values are in good agreement with neutron diffraction experiments (see Table 1). The dominant exchange coupling is J 3 Ni , which is much stronger than J 1 Ni . Note, that the experimental ratio of critical Mn . In both cases, J 2 is much smaller than the other two exchange couplings. In particular, J 2 and J 3 couplings might be considered as superexchange interactions involving the atoms in the path M−S1···S4−M for J 2 (see Figure 1d,e), where the S atoms are located in different sublayers, whereas the J 3 interaction is mediated by S atoms located in the same sublayer through the bridge M− S1···S2−M. One could expect stronger hybridization of the S p states and M 3d states within the same sublayer of S atoms. In addition, the calculations reveal that J 1 is AFM and FM for MnPS 3 and NiPS 3 , respectively. Note that for both MnPS 3   Effective Direct Exchange and AFM Couplings. In order to understand the origin of the exchange couplings and their relative strengths, we consider model Hamiltonians for the direct and superexchange interactions. Note, that the direct exchange discussed in this work is a (second order) kinetic exchange process, which involves only the hopping between the transition metal ions, with the ligands not explicitly involved, cf. ref 80, whereas a superexchange term is a fourth-(or higher-) order kinetic exchange process, which explicitly involves the hopping over the ligands, cf. ref 80. In addition, we are pointing out below the mechanisms which could impact the sign of J 1 .
To understand the origin of the magnetic couplings, we first write down a minimal many-body model which solely contains the valence electrons of the transition metal ion 81 the multiband Hubbard model. From this, using the second-order perturbation theory that is valid in the Mott insulating limit, we derive the (effective) direct exchange processes. By construction, all obtained spin couplings have to be AFM. Therefore, while surprisingly successful, the following simple analysis will not be able to explain the onset of the FM couplings in NiPS 3 (more on this at the end of the section).
In the monolayers of MnPS 3 and NiPS 3 , the metal atoms are surrounded by six sulfur atoms (MS 6 octahedron) and exhibit D 3d point group symmetry in trigonal anti-prismatic environment of the ligands (sulfur atoms), see Figure 2a. Hence, due to this trigonal crystal field effect, the d manifold splits in two disentangled subsets of bands. The only coupling between these subsets is the spin−orbit coupling. 82 The bands lower in energy are even (d x 2 −y 2 , d xy , and d z 2 ) with respect to the basal plane, while the bands higher in energy are odd (d xz and d yz ) with respect to the basal plane. The Mn ion is d 5 and Mn d bands split into half-filled even and half-filled odd bands; therefore, the magnetic coupling acquires contributions both from the even and the odd orbitals. Instead, because the Ni ion has a d 8 configuration, the even orbitals are fully occupied, and therefore, there is no magnetic contribution from these orbitals. Altogether, the minimal model is the five-band Hubbard model with a simplified structure of the Coulomb interactions (no spin on-site spin exchange and pair-hopping In this model, c imσ † creates an electron with spin σ = ↑, ↓ in a Wannier orbital |m⟩ = |x 2 − y 2 ⟩, |xy⟩, |xz⟩, |yz⟩, or |3z 2 − r 2 ⟩ at site i, and n imσ = c imσ † c imσ . ↑ (↓) indicates the spin up (down). The parameter t m,m′ i,i ′ is the hopping integral from orbital m at site i to orbital m′ at site i′. The on-site terms t m,m′ = ε m,m′ give the crystal field splitting. U and J H are the direct and (Hund) exchange terms of the screened on-site Coulomb interaction.  67 1.58 0.08 0.46 0.05 0.29 MnPS 3 (U = 4 eV, a = 6.00 Å, AFM-N) 21 0.4 0.03 0.15 0.08 0.38 MnPS 3 (U = 3 eV, a = 6.00 Å, AFM-N) 66 1.42 0.08 0.52 0.06 0.37 MnPS 3 (experiment, AFM-N) 64 1 were used. The AFM-N and AFM-z indicate the magnetic ground state of the system. The exchange couplings in the model calculations are obtained within the direct exchange mechanism [for MnPS 3 using eqs 2 and 4 whereas for NiPS 3 using eq 3]except for J 1 Ni for NiPS 3 , which contains also an important superexchange contribution and follows from eq 7; see text for further details. More details on the calculation of the magnetic exchanges within the DFT + U approach are present in Supporting Information. All values are given in meV.  where a and b are the dxz and d yz orbitals, respectively, while i and i′ are the Mn lattice sites. In this formula, we take in consideration both the cases with i = 1 ≠ i′ = 2 and i = i′ = 1. U Mn and J H Mn are the Coulomb repulsion and Hund coupling, respectively, in the case of the Mn atoms. The Hund's rule interaction between odd and even electrons yields a magnetic coupling between these electrons; therefore, the denominator depends on the occupancy of the even orbitals. By symmetry, the on-site energies are ϵ a

Mn2
. Similarly, we obtain the direct exchange constant for the valence electrons occupying the odd orbitals of the Ni atoms where a and b are the d xz and d yz orbitals, respectively, while j, j′ are the Ni lattice sites. In this formula, we take in consideration both the cases with j = 1 ≠ j′ = 2 and j = j′ = 1. U Ni and J H Ni are the Coulomb repulsion and Hund coupling, respectively, in the case of the Ni atoms. By symmetry the on- Ni2 . Finally, the direct exchange constant for the even Mn orbitals is as follows where c, d and e are the orbitals d z 2 , d x 2 −y 2 , d x,y respectively, while i and i′ are the metal lattice sites. In this formula, we take in consideration both the cases with i = 1 ≠ i′ = 2 and i = i′ = 1. U Mn and J H Mn are the Coulomb repulsion and Hund coupling, respectively, in the case of the Mn atoms. By symmetry the onsite energies are ϵ c

Mn2
. The obtained (total) direct exchange, J = J odd + J even , is positive for both Mn and Ni ion and for any distance j − j′. Hence, as already mentioned, it is always AFM by construction.
Numerical Evaluation of the AFM Couplings. The Wannier basis provides us with ab initio values of the hopping integrals and crystal field splitting. We calculate the hopping parameters and the on-site energies using the interpolated band structure of the Wannier functions of the d-subsector.  (3NN). In the case of odd orbitals, the 3NN couplings are greater than the 1NN couplings, even by an order of magnitude; therefore, it is very important to consider these hopping amplitudes in the calculations. On the other hand, for even orbitals, the 1NN couplings are greater than the 3NN couplings. The 2NN couplings are always smaller with respect to the 1NN couplings and 3NN couplings. In the case of even orbitals, we neglect the difference between the on-site energies assuming the following approximation ϵ c ). Now, we will numerically evaluate the second-and thirdneighbor direct exchange as a function of the first-neighbor direct exchange. For the odd subsector of MnPS 3 , we obtain J 2 Mn,odd = 0.037J 1 Mn,odd and J 3 Mn,odd = 2.026J 1 Mn,odd . For the even subsector of the MnPS 3 , we obtain J 2 Mn,even = 0.015J 1 Mn,even and J 3 Mn,even = 0.016J 1 Mn,even . Considering that J Ni,even = 0 due to fully occupied even orbitals, for the odd subsector of NiPS 3 , we obtain J 2 Ni,odd = 0.047J 1 Ni,odd and J 3 Ni,odd = 11.09J 1 Ni,odd . If we consider the total direct exchange value as the sum of the odd and even sector, we have  Table 1). Note that the values are overestimated in comparison to experimental values (see Table  1); however, the signs and the dominant contributions are the same.
When we numerically evaluate the direct exchange of NiPS 3 , we obtain J 1 Ni = (8.957 eV·meV)/(U Ni + J H Ni ) for the firstneighbor coupling, J 2 Ni = 0.047J 1 Ni and J 3 Ni = 11.09J 1 Ni for the second and third neighbors, respectively. Remarkably, the 3NN magnetic direct exchange is larger than the 1NN exchange in the odd case. Using a Coulomb repulsion of 6 eV for Ni and a Hund coupling of 0.7 eV, we obtain the numerical values equal to J 1 Ni = 1.3 meV, J 2 Ni = 0.06 meV, and J 3 Ni = 14.8 meV, where last two are reported in Table 1. Note, that J 2 Ni and J 3 Ni are in good agreement with experimental values (see Table 1), while J 1 Ni is significantly different and has a wrong sign. Therefore, in the case of J 1 Ni , a more complex superexchange model will be presented further in the text.
Altogether, we obtain that the simple direct exchange scheme gives J 3 Ni ≫ J 1 Ni ≫ J 2 Ni and J 1 Mn > J 3 Mn ≫ J 2 Mn . In the Ni case, the leading term is J 3 Ni , while in the Mn case, the leading term is J 1 Mn . The reason for this different behavior comes from the different filling that produces J Ni,even = 0. The calculated direct exchange couplings are qualitatively in agreement with the magnetic couplings obtained experimentally or using DFTexcept that J 1 Ni and J 2 Ni are FM due to more complex magnetic exchange not taken into account by the (simple) direct exchange scheme (see below). Even though the latter coupling is relatively small, note that considering in more detail, such discrepancy is relevant for an accurate description of the magnetic coupling and, hence, of the magnetic critical temperature. The different leading exchange terms in Mn and Ni compounds open the way to manipulate the magnetism by tuning the concentration of Mn and Ni compounds or by adding new magnetic materials as dopants.
Superexchange and FM Coupling in NiPS 3 . So far, we considered an effective direct exchange model, which solely contained the exchange processes due to the electrons hopping between the transition metal (Mn or Ni) ions. Note that such a model should be considered an effective one, for in reality, the The Journal of Physical Chemistry C pubs.acs.org/JPCC Article hopping between the neighboring Mn or Ni ions is predominantly mediated by the sulfur ions. Hence, a natural extension of the direct exchange model should explicitly contain the exchange processes on the sulfur ionsin fact, it is such a nearest neighbor superexchange model 84 that is studied below to explain the onset of the nearest neighbor FM exchange in NiPS 3 . Note that for consistency, we comment at the end of this subsubsection why the more complex superexchange model is not needed to understand the other, the firstand third-neighbor, a magnetic couplings in MPS 3  that is, the effective direct exchange model is enough. We introduce the nearest neighbor superexchange model 84 for NiPS 3 by considering two Ni ions connected via two sulfur atoms over two 90°bonds, see Figure 3. Note that, because Ni 2+ is in a d 8 configuration and sulfur ions are fully occupied, it is easier to consider the hole language as discussed below. Moreover, in what follows, we neglect the small trigonal distortions and we assume an octahedral crystal field with the division in t 2g and e g orbitals [with the coordinate system defined in such a way that the xy plane coincides with the plane formed by the nearest-neighbor transition metal ion and sulfur, see Figures 2b and 3]. Hence, we begin by considering a fully atomic limit without hopping (zeroth order of perturbation theory in the small kinetic energy) in which there are two holes localized in two distinct Ni e g orbitals (d x 2 −y 2 , d 3z 2 −r 2 ; in what follows assumed to be energetically degenerate, see also above) and two (lying higher by the charge transfer energy Δ) empty p orbitals (p x , py) on sulfur. As before, due to the strong Hund's rule J H , the two Ni 2+ holes form a high spin S = 1 state. Now, let us perform a perturbation theory in the kinetic energy (over Coulomb repulsion U and charge transfer energy Δ) and consider the possible exchange processeswhich are of two kinds: First, there are direct exchange processes between the nickel ions, see Figure 3a. By definition, these concern virtual occupancies of one of the nickel ions by three holes [with a relative energy cost of U + J H according to the simplified structure of the Coulomb interactions, see eq 1] and are possible once the hole can directly hop back and forth between the nickel orbitals under consideration. According to the Slater−Koster scheme, 85 which is qualitatively confirmed by our DFT calculations, the latter is allowed between the pair of d x 2 −y 2 orbital (ddπ hopping element) and between the pair of d 3z 2 −r 2 orbitals (ddσ/4 hopping element; the small ddδ element can be neglected), cf. Figure 3(a1,a2); note that the hopping between the d x 2 −y 2 and the d 3z 2 −r 2 orbital vanishes in this geometry, cf. Figure 3(a3). Altogether we obtain the nearestneighbor direct exchange contribution where we assumed a typical relation between the Slater− Koster hopping integrals ddσ ≈ 2ddπ. It is important to state here that the above direct exchange process is different than the ef fective direct exchanges process defined in the previous subsectionsfor the latter ones may include all indirect hoppings (e.g., via sulfur) between the two nearest neighbor nickel ions. Second, there are superexchange processes between the nickel ions, see Figure 3b. By definition these concern virtual Finite direct-AFM exchange processes due to (a1) nonzero hopping elements ∝ ddπ between the nearest-neighbor d x 2 −y 2 orbitals on nickel and (a2) nonzero hopping elements ∝ ddσ between the nearest-neighbor d 3z 2 −r 2 orbitals on nickel. Lack of direct-AFM exchange processes due to vanishing hopping elements between the nearest-neighbor d x 2 −y 2 and d 3z 2 −r 2 on nickel. Finite FM superexchange processes due to (b1) nonzero hopping elements ∝ pdσ between the nearest-neighbor d x 2 −y 2 orbitals on nickel and the nearest-neighbor p x /p y orbital on sulfur, (b2) nonzero hopping elements ∝ pdσ between the d 3z 2 −r 2 orbitals on nickel and the nearest-neighbor p x /p y orbital on sulfur, and (b3) nonzero hopping elements ∝ pdσ between the d x 2 −y 2 orbital on nickel and the nearest-neighbor p x /p y orbital on sulfur and the d 3z 2 −r 2 orbital on nickel and the nearest-neighbor p x /p y orbital on sulfur. See text for more details.
The Journal of Physical Chemistry C pubs.acs.org/JPCC Article occupancies of one of the sulfur ions by two holes (with an energy cost associated with the charge transfer 2Δ for antiparallel spins or 2Δ − J H for parallel spins b ) and are possible once the two holes hop back and forth between the sulfur and nickel orbitals under consideration. Again, using the Slater−Koster scheme, 85 which is negative (FM) due to the lowering of the energy for a FM (virtual) occupancy of sulfur by two holes in the superexchange process. Note that for a given Ni−Ni bond, there are always two superexchange processes: one over the top-left and one over the bottom-right sulfur (hence the factor of two within the square brackets in the above formula). Moreover, the process depicted in Figure 3(b3) has to be multiplied by a factor of two, for one can interchange the position of the d x 2 −y 2 and the d 3z 2 −r 2 orbital and in this way double the amplitude of this process. Finally note that, due to the fact that we have spins S = 1 on nickel, overall the above superexchange process is reduced by a factor 1/2 with respect to to an analogous one for the S = 1/2 on copper (the superexchange processes have to be projected on the high-spin S = 1 states on both nickel ionshence a factor of (1/ 2 ) 2 reduction).
Let us now comment why the above FM superexchange mechanism is not important for the nearest neighbor exchange in MnPS 3 . The reason for this is that in the case of manganese ions, we are a bit closer to the situation discussed in ref 84, which, for instance, shows the AFM exchange coupling in the case of half-filled t 2g subshells of Cr 3+ in LiCrS 2 . More precisely, the situation for the Mn 2+ ions is as follows. On one hand, the AFM direct exchange is much stronger for Mn 2+ because one of the t 2g electrons (the d xy ) can hop over the ddσ bond. On the other hand, the superexchange also contains an additional strong 84 AFM contribution due to the superexchange processes over one p z sulfur orbitalwhich strongly hybridizes with the two nearest t 2g orbitals. Altogether, as confirmed by the effective direct exchange studies in the previous subsection, these two mechanisms originating in the t 2g exchange processes easily overcome the (above-described) FM processes for the e g orbitals.
Finally, we mention that the (surprisingly) strong thirdneighbor coupling in MPS 3 can easily be explained using the superexchange model. First, for the t 2g sector, the strong AFM third-neighbor coupling is already discussed in ref 84see Figure 7 of the reference. Second, one can easily imagine that a similarly strong AFM coupling can also be realized for the e g electrons: in order to understand it, one just needs to replace the two third-neighbor d xy orbitals with the d x 2 −y 2 orbitals and rotate the sulfur p orbitals by 90°in the process shown in Figure 7 of ref 84. In fact, the latter one should have a very high amplitude and hence a really strong third-neighbor exchange in NiPS 3 .
Numerical Evaluation of the FM Coupling in NiPS 3 . Having derived the direct exchange and superexchange processes between the nearest neighbor nickel ions, we are now ready to estimate the contributions of these both (competing) spin interaction terms. To this end, we assume that Coulomb repulsion U = 6 eV (the chosen value of our DFT + U approach, see Table 1) and take a typical (for 3d transition metal compounds) value of Hund's exchange J H = 0.7 eV. Next, based on the DFT calculations, we estimate that (i) the charge transfer energy Δ = 3 eV, (ii) the hopping pdσ = 0.9 eV, and (iii) the hopping ddσ = 0.05 eV. From this and using eqs 5 and 6, we can easily calculate the AFM and FM contributions to the spin exchange As the second exchange is larger (by absolute value), we conclude that it is a relatively strong superexchange along the 90°nickel-sulfur-nickel bonds, which triggers the FM exchange along the nearest neighbor nickel ions. The sum of the aforementioned contributions leads to FM exchange coupling (J 1 Ni = −4.9), which is in good agreement with experimental studies (see Table 1).
At this point, we would like to comment on the fact that for the edge-sharing copper chains, typically one entirely neglects the direct exchange cf. ref 86. The reason for this is that, in the copper chains, the direct copper hopping ddσ = 0.08 eV, 87 whereas usually one assumes that for the cuprates, pdσ ≈ 1.5 eV. With such hopping parameters (and assuming the typical cuprate values of U = 8 eV, Δ = 3 eV, and oxygen J H = 0.7 eV), one can immediately see that for the copper chains with 90°g eometry, J Cu,direct ≈ 3 meV ≪ |J Cu,SE | ≈ 50 meVwhich justifies why the direct copper exchange is typically neglected in such studies. In reality, the FM exchange is actually really small in such copper chains due to the angle along the copper− oxygen−copper bond not being strictly equal to 90°.
A strong suppression of the direct exchange is also observed in CrI 3 , 88 which is an experimentally confirmed 2D ferromagnet. CrI 3 is isostructurally identical to MPS 3 compounds when P atoms are removed and has partially-filled t 2g shells, just as MnPS 3 . This might suggest that CrI 3 should predominantly have AFM interactions, just as MnPS 3 . However, this is not the case because the relative strength of the direct and superexchange mechanisms is different in CrI 3 and MnPS 3 . In CrI 3 , the Cr−I−Cr bond angle is very close to 90°, leading to a strong FM superexchange, while the Cr−Cr distance is relatively large (3.95 Å), giving rise to a very weak AFM direct exchange. Interestingly, due to the larger distances between the magnetic atoms in CrI 3 , the third-neighbor exchange is also less relevant.
Doped Systems (M 3/4 ,X 1/4 )PS 3 . Our strategy is to try to induce long-range FM order via chemical doping, keeping the two-dimensional structure of the mother compounds intact. In principle, the chemical doping could influence the magnetic exchanges, especially the long-range ones. To this end, we consider doping of the host systems with distinct possible 3d elements. In particular, at low doping concentration, the main exchange couplings are the magnetic exchanges between the host magnetic atoms (J i M for i = 1, 2, 3) and between the host The Journal of Physical Chemistry C pubs.acs.org/JPCC Article impurity atoms (J i MX for i = 1, 2, 3). Although for most of the 3d doping, the Mn−Mn and the Ni−Ni magnetic couplings remain AFM, in the case of Cr doping, the AFM exchanges between the magnetic atoms of the host turn into FM ones. 89 Thus, we performed DFT + U calculations for 25% concentration of the magnetic dopants (X = Cr, Mn, and Ni) in the hosts of NiPS 3 and MnPS 3 . We have employed various structural arrangements of the atoms (see Figure 4) and collinear spin configurations such as AFM Neel (AFM-N), zigzag (AFM-z), stripy (AFM-s), and FM case. For each of the configurations, the atomic positions have been fully optimized, keeping the lattice constants equal to the pure optimized monolayer structures (see Methods).
Let us first discuss the impact of the dopant on the energetic and structural properties of the hosts. The energy difference between a particular configuration and the magnetic ground state for each considered alloy is presented in Figure 5. Note that each of the employed alloys preserves the magnetic ground state of the host, independent of the structural arrangement of the impurity and the type of the dopants (see Figure 5). In addition, the favorite Mn and Ni dopant position is at 1NN (see Figure 5c,d), whereas the Cr dopants prefer to lie further apart, in particular, at 2NN (see Figure 5a) and 3NN (see Figure 5b) structural positions for NiPS 3 and MnPS 3 hosts, respectively. This reveals that the Mn and Ni dopants have tendency to cluster, while the Cr ions prefer to spread over the host. From now on, we only discuss the magnetic ground-state configurations of the alloys (the one which exhibits the lowest energy in Figure 6 for particular alloy). The magnetization and the band gaps of these systems are collected in Table 2. All considered alloys exhibit a semiconducting behavior (see Table 2). However, only (Ni 3/4 ,Cr 1/4 )PS 3 and (Ni 3/4 ,Mn 1/4 )PS 3 alloys have a nonzero net magnetization. Thus, we focus our discussion on these two ferrimagnetic systems. Owing to the different spins of the host and dopant atoms and particular arrangement of the dopants in the NiPS 3 host, the ferrimagnetic state appears in these systems. The band structure and the orbital projections of these two ferrimagnetic alloys are presented in Figure 6. Note that in the (Ni 3/4 ,Cr 1/4 )PS 3 case, bands close to the Fermi level are mainly composed of 3d states of Cr dopants, causing a sizeable reduction in the energy gap compared to pure NiPS 3 (2.3 eV for U = 6 eV). In the case of Mn impurities, the valence band maximum is mainly composed of Mn 3d states, whereas the conduction band minimum consists of very flat bands of Ni atoms (see Figure 6b).
Next, we examine the mixed exchange couplings between the metal atom of the host and the impurity J i MX within the classical Ising Hamiltonian in the honeycomb lattice for the smallest possible cell containing 25% dopant concentration (for details see Supporting Information and Figure S1 therein). Note that the FM exchange couplings are obtained for the Mn and Ni nearest neighbors and Cr−Cr ions at the secondnearest-neighbor distance (see J 2 X in Table 3). The similar results for Cr have been recently reported in ref 90. In addition, the different values for J 1 MnNi (−1 meV) and J 1 NiMn (−0.7 meV) stems from the different lattice parameters of the hosts (see Table 3). Although the mixed exchange coupling obtained for (Ni 3/4 ,Mn 1/4 )PS 3 correctly reflected the 1NN AFM-z ground-state configuration for this system, the difference between the 1NN AFM-N and 1NN AFM-z is just 1.5 meV per magnetic ion (see Figure 5c), thus resulting in magnetic frustration between the Mn and Ni atoms at the honeycomb sites due to the competition between the Neél and zigzag configurations, similar to the results reported in ref 35. Moreover, the critical temperature is directly related to the strength of the exchange couplings. Generally, the mixed exchange couplings are smaller than those in the corresponding metal atoms of the hosts. Thus, the critical temperature of the mixed structure is expected to be smaller than that for the corresponding pure system, which is in line with recent experimental reports on the Ni 1−x Mn x PS 3 alloy 35 and with series of the mixed-system studies, where the suppression of T N temperature with dopant substitution is reported. 27,29,34 In addition, one should expect the further reduction of T N temperature of the employed mixed systems due to the

■ CONCLUSIONS
First, we examined the magnetic and electronic properties of the AFM-ordered systems MnPS 3 and NiPS 3 (without magnetic impurities). We presented a qualitative explanation for the relative ratio of the different nearest neighbor (first, second, and third neighbor) exchange couplings by studying an effective direct exchange for both MnPS 3 and NiPS 3 . In particular, we demonstrated that the third-neighbor exchange dominates in NiPS 3 due to the filling of the t 2g subshell, whereas for MnPS 3 , the first-neighbor exchange is prevailed, owing to the presence of the t 2g magnetism. We showed in this work that the onset of the nearestneighbor FM coupling in NiPS 3 is due to the relatively strong (FM) superexchange, which is enabled by the complete filling of the t 2g shell and by the (close to) 90°nickel−sulfur−nickel nearest neighbor bond. Nevertheless, even these relatively "fortunate" circumstances do not guarantee that the strong further neighbor exchange can also be FM. The reason for this lies in the further neighbor bonds, which are no longer of a "90°variety"and the latter bonds are, apart from specific Jahn−Teller effects, the only way to obtain FM couplings in Mott insulators. From a more general perspective, this study confirms the paradigm that the AFM couplings are natural to the Mott insulating compounds.
Second, we examined the properties of the MnPS 3 and NiPS 3 compounds doped with 25% impurities (Ni in MnPS 3 , Mn in NiPS 3 , and Cr in both hosts). It turned out that all of the investigated alloys are Mott insulating, albeit with generally smaller band gaps than those of the corresponding host. Crucially, we demonstrated an extreme robustness of the AFM phases against impurity doping of the MnPS 3 and NiPS 3 compounds. Hence, ferromagnetism cannot be easily stabilized by impurity doping, in agreement with the above paradigm. Nevertheless, as the dopants have different spins than the pure phases, the alloys exhibit ferrimagnetic properties for particular arrangements of dopants.
The Mn and Ni impurities prefer to form dimers within the host, whereas the Cr dopants prefer to be further apart. Interestingly, unlike for the hosts, the first and second (dopant−host) exchange couplings are of similar order of magnitude. The latter leads to frustration in the case of AFM couplings. We suggest that this may be one of the reasons of the observed lower magnetic ordering temperature of the doped systems. 35 Our work sheds light on the origin of magnetism in the AFM family of transition metal phosphorus trichalcogenides by pointing out the mechanisms which govern the benchmark compounds, thus extending the fundamental knowledge of 2D magnetism.      Note that the finite Coulomb repulsion cost on sulfur just renormalizes the (anyway quite hard to estimate) charge transfer energy.